Critical behavior of self-assembled rigid rods on triangular and honeycomb lattices 
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Using Monte Carlo simulations and finite-size scaling analysis, the critical behavior of self- 
assembled rigid rods on triangular and honeycomb lattices at intermediate density has been studied. 
The system is composed of monomers with two attractive (sticky) poles that, by decreasing tem- 
perature or increasing density, polymerize reversibly into chains with three allowed directions and, 
at the same time, undergo a continuous isotropic-nematic IN transition. The determination of the 
critical exponents, along with the behavior of Binder cumulants, indicate that the IN transition 
belongs to the q — 1 Potts universality class. 
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I. INTRODUCTION 

Self-assembly has been considered for over 50 years to 
be central to understanding structure formation in living 
systems [l|. As a consequence, an increasing interest has 
been devoted to enhance our understanding of the the- 
oretical basis of the fundamental mechanisms governing 
self-assembly and the observables required to character- 
ize the interactions driving thermodynamic self-assembly 
transitions Despite these studies, the knowledge of 

how this process works is still incomplete and many of the 
basic principles characterizing this type of organization 
are unclear. 

It is obvious that a complete analysis of the self- 
assembly phenomenon is a quite difficult subject because 
of the complexity of the involved microscopic mecha- 
nisms. For this reason, the understanding of simple mod- 
els with increasing complexity might be a help and a 
guide to establish a general framework for the study of 
this kind of systems, and to stimulate the development 
of more sophisticated models which can be able to repro- 
duce concrete experimental situations. 

In this line of work, two previous articles [1, re- 
ferred to as papers I and II, respectively, were devoted 
to the study of a system of self- assembled rigid rods ad- 
sorbed on a two-dimensional lattice. In paper I, Tavares 
et al. studied a system composed of monomers with two 
attractive (sticky) poles that polymerize reversibly into 
polydisperse chains and, at the same time, undergo an 
isotropic-nematic (IN) continuous phase transition. So, 
the interplay between the self-assembly process and the 
nematic ordering is a distinctive characteristic of these 
systems. Using an approach in the spirit of the Zwanzig 
model [l(]| |. the authors found that nematic ordering en- 
hances bonding. In addition, the average rod length was 
described quantitatively in both phases, while the loca- 
tion of the ordering transition, which was found to be 
continuous, was predicted semiquantitatively by the the- 
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ory. With respect to the characteristics of the phase tran- 
sition, Tavares et al. assumed as working hypothesis that 
the nature of the IN transition remains unchanged with 
respect to the case of monodisperse rigid rods on square 
lattices, where the transition is in the 2D Ising universal- 
ity class [Ill4l3| . 



Paper II was a step further, analyzing the universal- 
ity class of the IN phase transition at intermediate den- 
sity. For this purpose, an extensive work of Monte Carlo 
(MC) simulations and finite-size scaling (FSS) analysis 
was carried out. The obtained results showed that the 
self-assembly process affects the nature of the transition. 
Thus, the accurate determination of the critical expo- 
nents indicated that, for a square lattice, the universality 
class of the IN transition changes from 2D Ising-type for 
monodisperse rods without self-assembly [HI to q = 1 
Potts-type for self-assembled rods. 



In this context, the objectives of the present paper are 
(1) to extend the previous work to triangular and hon- 
eycomb lattices using the same techniques developed in 
paper II; and (2) to study the effect of the lattice struc- 
ture on the critical behavior of self-assembled rigid rods. 
For this purpose, MC simulations [3] supplemented by 
analysis using FSS theory [l5| have been carried out to 
study the critical behavior in a system of self-assembled 
rigid rods deposited on triangular and honeycomb lat- 
tices. As in paper II, the calculations were developed 
at constant temperature and different densities, being 
(8/9 c — 1) the normalized scaling variable, where 9 and 6 C 
represent density and critical density, respectively. The 
accurate determination of the critical exponents, along 
with the behavior of Binder cumulants, confirmed that 
the IN transition of self-assembled rigid rods on triangu- 
lar and honeycomb lattices at intermediate density be- 
longs to the q = 1 Potts universality class. The outline 
of the paper is as follows. In Sec. II we describe the 
lattice-gas model and the simulation scheme. In Sec. Ill 
we present the MC results and the general conclusions. 
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FIG. 1: (a) Schematic representation of the set of vectors {fti, H2, n^} for triangular and honeycomb lattices, (b) Self-assembled 
rigid rods on triangular lattices, (c) Self-assembled rigid rods on honeycomb lattices. See discussion in the text, (d) Available 
configurations for a monomer on a honeycomb lattice. 



II. LATTICE-GAS MODEL AND MONTE 
CARLO SIMULATION SCHEME 

As in papers I and II, we consider a system of self- 
assembled rods with a discrete number of orientations 
in two dimensions. In this case, the surface is repre- 
sented as an array of M = L x L adsorption sites in 
a triangular or in a honeycomb lattice arrangement and 
periodic boundary conditions. N particles are adsorbed 
on the substrate with three possible orientations along 
the principal axis of the lattice. These particles interact 
with nearest-neighbors (NN) through anisotropic attrac- 
tive interactions. Thus, a cluster or uninterrupted se- 
quence of bonded particles is a self-assembled linear rod. 

Under these considerations, the adsorbed phase is char- 
acterized by the Hamiltonian 

H = -w ^ ' °i ■ °j I div 1, (1) 

where indicates a sum over NN sites; div means 

integer division; w represents the NN lateral interaction 
between two neighboring i and j, which are aligned with 
each other and with the intermolecular vector r^ ; and 
<7i is the occupation vector with (Tj = if the site i is 
empty, = £\ if the site i is occupied by a particle 
with orientation along the cci-axis, Bi = £2 if the site i is 



occupied by a particle with orientation along the X2-axis, 
and di = £3 if the site i is occupied by a particle with 
orientation along the X3-axis. 

At fixed temperature, the average rod length increases 
as the density increases and the polydisperse rods will 
undergo an ncmatic ordering transition [8, 9]. In order 
to follow the formation of the nematic phase, we use the 
order parameter defined in previous work 0, Q , which in 
this case can be written as 

s = \ni+n 2 + ri 3 \ 
\ni\ + \n 2 \ + |n 3 | 

where each vector n m is associated to one of the 3 possible 
orientations (or directions) for a chain on the lattice. In 
addition, (1) the Hi's lie in the same plane (or are co- 
planar) and point radially outward from a given point 
P which is defined as coordinate origin; (2) the angle 
between two consecutive vectors, rJj and n i+ i, is equal 
to 27r/3; and (3) the magnitude of n, is equal to the 
number of fc-mers aligned along the i-direction [see Fig. 
1(a)]. Note that the Hi's have the same directions as the q 
vectors in Ref. fl6l |. These directions are not coincident 
with the allowed directions for the chains on the real 
lattice. 

The concept of linear rod is trivial for triangular lat- 
tices [see Fig. 1(b)]. However, in a honeycomb lattice, 
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the geometry does not allow for the existence of a linear 
array of monomers. In this case, we call "linear rod" to 
a chain of adjacent monomers that can be assembled in 
only three types of sequences, defining three directions in 
similar way to the triangular lattice [see Fig. 1(c)]. Two 
"types" of sites can be recognized in the honeycomb lat- 
tice: A (grey) and B (black). Monomers (colored accord- 
ing to the adsorption site) can be placed on these sites 
in three different orientations: Ax, A%, and A3 (for the 
A-sites); and B\, B 2 , and B 3 (for the B-sites) as shown 
in Fig. 1(d). So, the self-assembled chains are formed 
by alternating A and B monomers. In the particular 
case shown in Fig. 1(c), C\, C2, and C3 are built by the 
following sequences: B\A\B\A\B\A\ (Ci), B2A2B2A2 
(C2), and A3B3A3B3A3 (C3). The figure also shows an 
unassociated B-monomer (-B2). 

The problem has been studied by canonical Monte 
Carlo simulations using an vacancy-particle-exchange 
Kawasaki dynamics [l7[ and Metropolis acceptance prob- 
ability IS] . Typically, the equilibrium state can be well 
reproduced after discarding the first 10 7 Monte Carlo 
steps MCS. Then, the next 10 9 MCS are used to com- 
pute averages. All calculations were carried out using 
the parallel cluster BACO located at Instituto de Fi'sica 
Aplicada, Universidad Nacional de San Luis-CONICET, 
San Luis, Argentina. This facility consists of 60 personal 
computers each with a 3.0 GHz Pentium-4 processor and 
90 personal computers each with a 2.4 GHz Core 2 Quad 
processor. 

In our Monte Carlo simulations, we set the tempera- 
ture T, varied the density 9 = N/M and monitored the 
order parameter (5, which can be calculated as simple av- 
erage. The quantities related with the order parameter, 
such as the susceptibility x, and the reduced fourth-order 
cumulant Ul introduced by Binder fl4j | were calculated 
as: 
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FIG. 2: Size dependence of the order parameter as a function 
of coverage. Upper-left inset: Order parameter in the regime 
of high densities. Lower-right inset: Size dependence of the 
susceptibility as a function of coverage. The error in each 
measurement is smaller than the size of the symbols, (a) 
Triangular lattices, (b) Honeycomb lattices. 
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U L = 1 
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where (■ ■ • ) means the average over the MC simulation 
runs. 

In addition, in order to discuss the nature of the phase 
transition, the fourth-order energy cumulant Ue was ob- 
tained as 



U E = 1 



3{# 2 ) 5 



(5) 



III. RESULTS AND CONCLUSIONS 

As discussed in section I, the phase diagram of a sys- 
tem of self- assembled rigid rods on square lattices has 



been recently reported by Tavares et al. |8j. The au- 
thors showed that the critical density, at which the IN 
transition occurs, increases monotonically as ksT/w is 
increased. Thus, the nematic phase is stable at low tem- 
peratures and high densities [see Fig. 1(a) in Ref. 0]. 
Later, and based on this finding, the nature of the IN 
transition at intermediate densities was studied in Ref. 
Q. There, the calculations were developed at w = iksT 
and different densities. 

In the present study, we set the lateral interaction to 
w = 4.5/csT [a close value to that used in Ref. @], 
thus allowing a direct comparison with previous results 
for square lattices. Accordingly, the density was varied 
around half coverage. For each value of 6, the effect of 
finite size was investigated by examining lattices with L 
ranging from 30 to 180. 

We start with the calculation of the order parameter 
(Fig. 2), and susceptibility (lower-right insets in Fig. 
2) plotted versus 8 for several lattice sizes as indicated. 
Curves in part (a) [(b)] correspond to triangular [hon- 
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FIG. 3: Curves of Ul(9) vs # for lattices of different sizes. 
From their intersections one obtained 6 C . In the lower-right 
inset, the data are plotted over a wider range of densities. 
Upper-left inset: Coverage variation of Ue for various lattice 
sizes, (a) Triangular lattices, (b) Honeycomb lattices. 



eycomb] lattices. As it can be observed, 8 appears as 
a proper order parameter to elucidate the phase tran- 
sition. When the system is disordered, all orientations 
are equivalents and 6 is zero. In the critical regime, the 
particles align along one direction and 6 increases con- 
tinuously to one, remaining constant up to full coverage 
(see upper-left insets in Fig. 2, where the order parame- 
ter is shown up to 9 = 1). In other words, nematic order 
survives until 9 = 1. This finding allows us to discard 
the existence of a reentrant nematic transition at high 
densities as speculated in Ref. Q. With respect to the 
susceptibility, the curves show a single peak which grows 
and sharpens as the lattice size is increased. 

The critical density has been estimated from the plots 
of the reduced fourth-order cumulants Ul(9) plotted vs. 
9 for several lattice sizes (see Fig. 3). In the vicinity of 
the critical point, cumulants show a strong dependence 
on the system size. However, at the critical point the cu- 
mulants adopt a nontrivial value U*; irrespective of sys- 
tem sizes in the scaling limit. Thus, plotting Ul(9) for 
different linear dimensions yields an intersection point 
U* , which gives an accurate estimation of the critical 



density in the infinite system. In this case, the values 
obtained for the critical density were 6 C = 0.503(7) for 
triangular lattices [part (a)] and 9 C = 0.659(6) for hon- 
eycomb lattices [part (b)]. 

With respect to the value of the cumulant at the tran- 
sition density, which allows us to make a "preliminary" 
identification of the universality class of the transition 
HI EHHl, we obtained U* = 0.639(2) for triangular 
lattices [part (a)] and U* = 0.637(7) for honeycomb 
lattices [part (b)]. These values of U* are practically 
indistinguishable from previous estimates for the two- 
dimensional one-state Potts model (see for instance Ref. 
H, where U* = 0.639(3)). This result may be taken as 
a first indication of universality. In the lower-right inset, 
the data are plotted over a wider range of densities. As 
can be seen, the curves exhibit the typical behavior of Ul 
in the presence of a continuous phase transition. Namely, 
the order-parameter cumulant shows a smooth increase 
up to its ordered-phase value of 2/3 instead of the char- 
acteristic deep (negative) minimum, as in a first-order 
phase transition [ljj . 

In order to discard the possibility that the phase tran- 
sition is a first-order one, the energy cumulants [Eq. (9)] 
have been measured. As is well known, the finite-size 
analysis of Ue is a simple and direct way to determine 
the order of a phase transition [2l2H2~4j ] . Upper-left inset 
in Fig. 3 illustrates the energy cumulants plotted versus 
9 for different lattice sizes. Data in part (a) correspond to 
triangular lattices. Data in part (b) correspond to honey- 
comb lattices. As is observed, Ue has the characteristic 
behavior of a continuous phase transition. Namely, the 
curves of Ue show a dip close to the critical density for 
all system sizes, but this minimum scales to 2/3 in the 
thermodynamic limit. This indicates that the latent heat 
is zero in the thermodynamic limit, which reinforces the 
arguments given in the paragraphs above. 

The results obtained up to here allow us to confirm the 
existence of a continuous phase transition at intermediate 
coverage and low temperature. In addition, the evalua- 
tion of the fixed point value of the cumulants indicates 
that, as in the case of square lattices [9|, the observed 
phase transition belongs to the q — 1 Potts universality 
class. To corroborate these findings, the critical behavior 
of the present model has been investigated by means of 
finite-size scaling (FSS) analysis. The FSS theory implies 
the following behavior of 6, x an( i Ul at criticality: 

b = L- p, l v l{L l l v t) 1 (6) 

X = L^» X (L 1/lJ e) (7) 

and 

U L - U L {L 1,u e), (8) 

for L — > oo, e — > such that L x l v (~ finite, where 
e = 9/9 c — 1. Here j3, 7 and v are the standard critical 
exponents of the order parameter (6 ~ — e* 3 for e — > 0~, 
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FIG. 4: (a) Data collapsing of the order parameter, SL 13 ^ vs 
lelL 1 /", and of the cumulant, Ul vs eL 1 ^ (inset) for triangu- 
lar lattices. The plots were made using 9 C = 0.503 and the 
exact percolation exponents v = 4/3 and j3 — 5/36. (b) Same 
as part (a) for honeycomb lattices and 6 C = 0.659. 
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FIG. 5: (a) Data collapsing of the the susceptibility, \L 
vs eL 1 ^ , for the curves in the inset of Fig. 2 (a). The plot was 
made using 8 C — 0.503 and the exact percolation exponents 
v = 4/3 and 7 = 43/18. (b) Same as part (a) for the curves 
in the inset of Fig. 2 (b) and 6 C = 0.659. 



L — > 00), susceptibility(x ~ |e| 7 for e — > 0, L — > 00) 
and correlation length £ (£ ~ |e|~^ for e — » 0,1/ — > 00), 
respectively. <5, \ an d are scaling functions for the 
respective quantities. 

According to eqs. ([HUH]), the scaling behavior can be 
tested by plotting {5)l/l v vs \e\L^ u , xL^l v vs eL 1 ^, 
and Ul vs eL 1 /" and looking for data collapsing. Using 
the exact values of the critical exponents of the ordinary 
percolation (one-state Potts model) v = 4/3, j3 = 5/36, 
and 7 = 43/18, we obtain an excellent scaling collapse as 
it is shown in Figs. 4 and 5. Data in part (a) correspond 
to triangular lattices and 6 C — 0.503. Data in part (b) 
correspond to honeycomb lattices and 6 C — 0.659. The 
study in Figs. 4 and 5 corroborates that the transition 
observed corresponds to the universality class of the q = 1 
Potts model. 

Finally, it is important to remark that the universality 
observed for self-assembled rigid rods on triangular and 
honeycomb lattices coincides with that reported in Ref. 
@ for square lattices. In other words, the universality 
class of the present model does not depend on the lattice 
geometry. This result (1) reinforces the the arguments 



given in Ref. |9j linking the thermal phase transition (IN 
phase transition) occurring in the system and the perco- 
lation behavior of the clusters of aligned monomers [25j j ; 
and (2) reveals a significant difference between the model 
here discussed and that of monodisperse rigid rods ad- 
sorbed on two-dimensional lattices with a discrete num- 
ber of orientations [TJ Ell HIL where the universality 
class of the IN phase transition occurring at intermedi- 
ate densities belongs to the 2D Ising universality class for 
square lattices and to the three-state Potts universality 
class for triangular and honeycomb lattices. 

In summary, we have used Monte Carlo simulations 
and finite-size scaling theory to study the critical prop- 
erties of self-assembled rigid rods on triangular and hon- 
eycomb lattices at intermediate density. The existence of 
a IN continuous phase transition was confirmed. In addi- 
tion, the scaling behavior of the system revealed that, as 
in the case of square lattices, the phase transition belongs 
to the q = 1 Potts universality class. 
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